Stochastic Cluster Series expansion for quantum spin systems 
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In this paper we develop a cluster- variant of the Stochastic Series expansion method (SCSE). 
For certain systems with longer-range interactions the SCSE is considerably more efficient than the 
standard implementation of the Stochastic Series Expansion (SSE), at low temperatures. As an 
application of this method we calculated the T — 0-conductance for a linear chain with a (diagonal) 
next nearest neighbor interaction. 
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Introduction — The development of efficient Quantum 
Monte Carlo (QMC) methods with loop-updates, like 
the original loop-algorithmi and the SSES^*^, has been 
a major advancement. They are very efficient for the 
anistropic Heisenberg models, hke the xxz-cha.m and can 
be generalized to more complicated Hamiltonians, but in 
some cases only with reduced performance. 

Here we study the xxz-cham with a diagonal next near- 
est neighbor interaction. This model is better suited 
for the description of fermionic systems since it takes 
into account that the interaction (Coulomb-repulsion) 
is long ranged. Furthermore, it is one of the simplest 
non-integrable systems. Consequently, this model has 
attracted the attention of many authorsjSi^ 

We find that a standard SSE-implementation performs 
only poorly on our model system. The reason for this lies 
in the fact that a certain transition in the operator loop 
update — called "bounce" — is given a relatively large 
weight. Syljuasen and Sandvik have shown that such a 
large bounce may affect the efficiency of the algorithm — 
especially at low temperatures. Hence, one should strive 
to find a way to reduce the bounce weight. 

Here we report, that a new implementation of the SSE- 
algorithm, the cluster-SSE, yields a considerably smaller 
bounce weight than the standard SSE, improving thus 
the efficiency considerably at low temperatures. 

Conventional SSE — We will now briefly discuss the 
conventional SSE-implementation and point out its diffi- 
culties. Our model of interest is a frustrated chain with a 
diagonal next-nearest neighbor interaction: (see Fig.^ 
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Note that only the interaction part is frustrated such 
that a Monte Carlo simulation will not suffer from the 
sign problem. Following Ref. 2 we start by splitting 
the Hamiltonian into a set of local operators, i.e., H = 
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where the constants C and C2 are needed to ensure that 
all matrix elements between S'^-eigenbasis-states are pos- 
itive. Using the Taylor expansion, the partition function 
may be written as 
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where M e N, : {1,...,M} —>■ f), and the am run 
over all 5'^-eigenbasis states with periodic boundary- 
conditions \ai) |aA/+i) along the imaginary-time axis. 
The factors (?!)(m)Q!m_(.i) in Eq. are called pla- 
quettes and are non- negative due to Hence, we can 
obtain the partition function by sampling the terms in 
Eq. with their relative weight factors over all spin- 
configurations \am) and local operators (j>{m) e t). 

Since each operator (f)(rn) acts only on two sites we 
may write 

l^"i,"2(s«i'«"2''Sni :S«2 ) := (am,0(rn)am+i) 

where s™ = S^\am) are spin variables and Wni,n2 is 
given by the prefactors of l|2Jl. Note that because of 
translation invariance Wni,n2 depends actually only on 
712 — ni . For brevity we will in the sequel denote the four 
arguments, which we call plaquette legs, by one super- 
argument, which we call the plaquette state. The SSE 
knows two updates: 

• In the diagonal update an insertion/removal of one 
plaquette is considered. 

• The loop update (constructs and) flips a sub-set of 
the spin variables such that the new conflguration 
is allowed — that means has non-zero weight. 

The loop update may be achieved step by stepi^ In each 
step we consider only one plaquette, whose state will be 
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changed from i to j by flipping two legs h and Ij. The 
"in-going" leg li is given and we have to choose the "out- 
going" leg Ij with a certain probability p{i — > 
Of course, the probability for a transition i j must be 
zero, if one of the states i or j has zero weight. 

These probabilities — yet to be determined — will play 
a central role in what follows. Therefore, we will explain 
in detail how they are obtained. 

They may equivalently be represented by (2^ • 4) x (2* • 
4)-matrices a"!'"'^ with entries 

The fact that the p{i — > j, U, Ij) are probabilities imposes 
for each li and i the following conditions: 

^^n.«.(*)=E«;:r)U)- (4) 

The detailed balance condition is satisfied if the matrices 
Q^ni,n2 g^j-g symmetric. 

Closer inspection shows that the matrices a"^'"^ are 
block diagonal where the dimension of the blocks d equals 
half the number of allowed plaquette states. Hence, if 
712 — ni = 1, we have d = 3 (as for the a;a;z-chain^; it 
would be four if pair creation and pair annihilation were 
not forbidden'^) whereas d = 2 for 712 — ni = 2. 

All the elements in one single block of a"i'"2 (differ in 
at least one state index. Hence, we may drop the leg 
indices h and Ij without causing ambiguities. Moreover, 
we will also omit the explicit indication of the site indices 
TT-i and 712 in a and W. 

For the blocks of length 3 a formal solution for the non- 
diagonal entries which meets detailed balance and 
Eqs. Q can easily be stated^ 

a,, = [W{i) + WU) - W{k)]/2 + [akk - a,, - a,,]/2 (5) 

where k is the third plaquette state in the same block of 
a as i and j. One sees that the diagonal entries of the 
matrices a remain free — apart from the final restriction 
that allprobabilities (entries of a) need be positive. 

Ref. (9 tells us how to dispose of the remaining degrees 
of freedom: The diagonal entries of a (baptized bounces^ 
since they correspond to the choice Ij = k) are always 
allowed, but turn out to be inefficient as they impede the 
update scheme. As a general ruley^ one can say that the 
most favorable among all possible solutions to Eqs. Q is 
the one with minimal diagonal entries. 

In the blocks of length three all bounce weights may 
be put to zero for a wide parameter range, ^ but in the 
blocks of dimension two Eqs. dictate one of the two 
bounces an and ajj to be equal to |VF(z) — M^(j')| = Jz2/2. 
A situation which is far from optimal. 

Cluster variant (SC'SE) — The problem with the SSE 
method is that some plaquettes — associated with opera- 
tors /i'''' — have only four allowed (diagonal) states. This 
can be avoided by splitting the Hamiltonian into small 
clusters instead of mere two-sites operators. We now in- 
troduce the Stochastic CZitster Series expansion (SCSE). 



13 5 7 




FIG. 1: The model Hamiltonian that will be discussed in this 
paper. Solid lines indicate a full Heisenberg-like interaction 
between the sites; dashed lines stand for sites coupled only by 
a z-z-term (Ising-like interaction). Two clusters — as defined 
in the text — are indicated by thicker lines. 
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FIG. 2: Pictorial representation of the plaquettes similar to 
the one used in Ref. ^ The operator itself is indicated by 
the horizontal central bar. The symbols stand for the six 
legs of the plaquettes. Equal symbols mean equal spins, and 
filled and empty symbols mean opposite spins. In this way 
the left diagram represents the eight diagonal plaquettes, and 
the right diagram, four non-diagonal plaquettes. The remain- 
ing non-diagonal plaquettes may be obtained by reflection 
with respect to the vertical central plaquette-axis and are not 
shown here. 

In the case of the frustrated chain we split the Hamilto- 
nian not into two- sites but three-sites operators (see 
Fig.©: 

^n,ri+l,ri+2 ~ •^aJ'^n 

^n,n+l.n+2 ~ JxS^ '5'„-|-i/4 

^n,n+l.n+2 = •-'2;>S'„^i5'„_|_2/4 

^n,n+l,n+2 ^ '^2;'^n-|-l'^n 4-2/4 

^n,n+l,n+2 = ^ ^- I'^i^n^n+l ^ ^n+l^n+2)~^ 
+ Jz2SnSn+2 + ^■ 

Note the factor 1/4 instead of 1/2 for with t = 
1,2,3,4. It stems from the fact that each non-diagonal 
operator is "distributed" between two clusters. The pos- 
sible plaquette states are shown in Fig. [3 

The diagonal update and loop construction remain un- 
changed with respect to the SSE. In the loop update the 
only difference is that W is now a function of 6 variables. 
Let us now turn to the matrix a. The transitions of the 
form h^''\ when t £ {1,2} and r e {3,4} will 

be ruled out from the beginning (for simplicity and be- 
cause we do not expect that they will give assistance in 
minimizing the bounce.) All other transitions may have 
positive probability. 

For a given plaquette — associated with an operator 
and a given in-going leg, whose site index is 77^ e 
{rj, 77 -|- 1,77 -|- 2}, we can directly tell the dimension of 
the corresponding block of a and give a solution for its 
entries. There are three cases to be considered: 

I) For t — 1,2 and = 77 -I- 2 (or i = 3, 4 and 77^ — n) 
we get only two equations from Eqs. (0J. Since both pla- 
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FIG. 3: Example of a block of dimension 2 (top) and 3 
(bottom) of Eqs. (01 • (Notation see Fig. |21) In- and possible 
out-going legs are indicated by arrows and connected by a 
line. (Note that in-going legs may always be out-going legs.) 
Upon flipping the two connected legs and interchanging in- 
and out-going leg one diagram becomes another. The other 
blocks of the same dimension may be obtained by reflection 
operations of the plaquettes. 
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FIG. 5; We compare simulations (Jz — 2J2z, 192 sites, 
T = OMJ^/ks, 2-10'^ MC-sweeps) for SSE and SCSE. Shown 
is the ratio of the time consume of the SSE and SCSE sim- 
ulations as well as the ratio of the mean (over the first 10 
Matsubara frequencies) statistical error of SSE and SCSE. 
X — N/2 and y are defined in Eq. Q. 
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FIG. 4: Two examples of a block of dimension 4 of Eqs. Q 
similar to Fig. 13 The other blocks of the same dimension 
may be obtained by reflection with respect to the horizontal 
or vertical axis. 

quette states (corresponding to non-diagonal /i*^*') have 
equal weights, we can put the bounce weight to zero. The 
proceeding of the loop is then deterministic, (see top of 
Fig. El) 

II) Otherwise if (i = 5, n.i ^ n+1) or {t ^ 5, rii ~ n+1) 
we have to consider three equations which may be treated 
as in the SSE. (see Eqs. (O and bottom of Fig. El) 

III) In the remaining cases we need to solve a sub- 
system of Eqs. of not less than four equations, (see 
Fig.^) Among the four plaquette states involved in this 
block we find two states corresponding to non-diagonal 

— which we call in and j„. From Fig.^lwe infer that if 
in corresponds to t G {1,2}, jn corresponds to t £ {3, 4}, 
et vice versa. The remaining two states correspond to 
diagonal operators; we call them id and jd- The entries 
of a are given by 

a,„,,, = ma^{[Wiid)/2 + Wiin) -W{jd)/2]/2,0} 
a,„,,, = mm{[Wijd)/2 + Wiin)-Wiid)/2]/2,W{in)} 
a,,,,, = min{[W{td) + W{jd)-2W{in)]/2,Wiid)} 
a,,,,, = max{0,Wijd)~Wiid)~2W{in)} 
aj„,fe = ai,.,fc, k = id,3d- 



Note that we assumed W{id) < W{jd) and exploited 
W{in) = W{jn)- With this choice only one diagonal 
entry (namely, cij^j^) may be non- negative. To be con- 
crete: only the situation depicted in the upper part of 
Fig.Hadmits a^^j^ = (J^ - J^)/2 > if > J^;. 

The SCSE is more intricate (one has to consider sev- 
eral cases separately) and the bounce cannot always be 
avoided, but there is a considerable improvement (at low 
temperatures) with respect to the SSE. 

Numerical results — In this paper we use the method 
proposed in Ref. 9 to calculate the conductance for the 
Hamiltonian We evaluate the expression 

giuJAl) ^ -UJAl/fl COs{ujMT){PxPy{iT))dT (6) 

Jo 

at the Matsubara frequencies ujm — 27TM{(3h)~^ , M e 
N. The conductance is obtained by extrapolating g{ijj) 
from the Matsubara frequencies to w = 0. The op- 
erators in Eq. ^ are given by (e is the charge unit) 
Px = eX]„>a; ^n- The uj = 0-value of g does not depend 
on X or 

To illustrate the improvement gained by introducing 
SCSE we performed for a special choice of parameters 
{Jz = 2J22) simulations oi g{ujM)- In Fig. Elwe plotted 
the ratios of the time consume along with the ratio of 
the statistical error in g{ojM) (for Af = 1, . . . , 10). The 
SSE is slightly faster (about 80 percent) but while the 
statistical error grows linearly with for SSE, it grows 
more modestly for the SCSE. 

Frustrated system — Using the Lanczos-method Zhu- 
ravlev, Katsnelson and coworkers^ obtained a very com- 
plete picture of the frustrated system given by the 
Hamiltonian (^). They set-up a low-temperature phase 
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diagram with two gapped and a gap-less phase. In 
the gap-less phase the system is — for a large range of 
parameters — very well-described by the Luttinger liquid 
picture.* In Fig.|nig(w = 0) — extrapolated by a quadratic 
fit — is plotted versus Jz2 for various Jz < Jx- For this 
region in parameter space Ref. finds a phase boundary 
between the gap-less and the gapped phase at Jz2 = Jx- 

In a Luttinger liquid the conductance equals the Lut- 
tinger parameteriSiiSiii and may hence be viewed as an ef- 
fective interaction. Since Jz2 mediates a nearest neighbor 
attraction and thereby reduces the effective interaction, 
the conductance first grows with Jz2 and then assumes 
its maximum approximately when 2Jz2 = Jz is satisfied. 
We see that within error bars the conductance does not 
vary on the phase boundary. 

As we work at a finite temperature of ksT — O.OlJx, 
the conductance goes smoothly to zero in the gapped 
phase such that a determination of the phase boundary 
from Fig. El becomes difficult. 
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FIG. 6: Conductance of the Heisenberg chain (400 sites, T = 
O.OlJx /ks) with next nearest neighbor interaction Jz2- The 
phase boundary from Ref. ^ is displayed. (We use OBC's, 
2 • 10^ MC-sweeps.) 

A remarkable fact of the phase diagram in Ref. '3 is the 
fact that the gap-less phase is not bounded in parame- 
ter space; it contains, e.g., the line with 2Jz2 = Jz- It 
is therefore interesting to study the behavior of the con- 
ductance on this line. However, if we increase Jz and Jz2 
simulating the conductance becomes more difficult. (The 
statistical error grows.) 

Only the SCSE allows us to compute the conductance 
for as large interaction values as Jz = 5, Jz2 — 2.5. The 
conductance is plotted in Fig.[7| Beside the statistical er- 
ror we have an extrapolation error depending on our ex- 
trapolation scheme. In Fig. [T] we compared two schemes: 



a quadratic fit to the first three Matsubara frequencies 
and a linear fit from the first six Matsubara frequencies. 
The former should give a smaller extrapolation error, but 
it enhances the statistical error of g{uj)- The latter has 
a larger extrapolation error, but it suppresses the sta- 
tistical error of 17(0;). For small parameters the two fits 
almost coincide. But for larger parameter values — when 
the statistical error increases — the quadratic fit starts to 
fluctuate. 

For this system the DC-Drude weight D and the 
susceptibility % at T = may be obtained by ex- 
act diagonalization^ The same is true for the conduc- 
tance via the relation g = ir^/Dx valid for a Luttinger 
liquidi-2ii^ Fig. [7| provides also the exact diagonalization 
data showing good agreement with the SCSE data within 
error bars. 

Conclusion — In this letter we showed that the per- 
formance of the SSE can be substantially improved by 
selecting a different splitting of the Hamiltonian. We 
used this strategy for an a;a::z-chain with next-nearest- 
neighbor-interaction, but it should also apply to systems 
which have additional plaquettes with a reduced num- 
ber of plaquette states. These systems include fermionic 
chains with arbitrary interaction part, which may also be 
coupled by an interaction term as well as (not necessar- 
ily one-dimensional) spin systems that have at the same 
time Heisenberg-like and Ising-like interaction bonds. 
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FIG. 7: Conductance of the Heisenberg chain (400 sites, 
T = O.OlJx /kB) with next nearest neighbor interaction Jz2 
along the line in parameter space where 2/22 = Jz- Cir- 
cles: quadratic extrapolation from the first three Matsubara 
frequencies. Squares: linear extrapolation from the first six 
Matsubara frequencies. (We use OBC's, 2 • lO'^ MC -sweeps.; 
For comparison Exact diagonalization results are given. 
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